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Abstract 

Virtual quark pair screening leads to breaking of the string between funda- 
mental representation quarks in QCD. For unquenched four dimensional lattice 
QCD, this (so far elusive) phenomenon is studied using the recently developed 
truncated determinant algorithm (TDA). The dynamical configurations were 
generated on an Athlon 650 MHz PC. Quark eigenmodes up to 420 MeV are 
included exactly in these TDA studies performed at low quark mass on large 
coarse (but O(a^) improved) lattices. A study of Wilson line correlators in 
Coulomb gauge extracted from an ensemble of 1000 two-flavor dynamical con- 
figurations reveals evidence for flattening of the string tension at distances 
R> 1 fm. 
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1 Introduction 



The behavior of the static energy of a quark-antiquark pair at large distance provides perhaps 
the most striking quahtative difference between quenched and fuU QCD. In the quenched 
theory virtual quark-antiquark pairs unconnected to external sources are neglected and 
the string tension rises indefinitely at large distance, whereas the full theory automati- 
cally screens the quark-antiquark potential at large distances by populating the vacuum 
with dynamical quark pairs. It is therefore not surprising that the demonstration of string 
breaking has long been regarded as a classic bellwether for testing the efficacy of dynam- 
ical QCD algorithms. Unfortunately, despite numerous studies 0| and the expenditure of 
a large amount of computational effort, the direct observation of string breaking in zero- 
temperature 4-dimensional unquenched QCD, in which the string tension is seen to become 
essentially flat at large distances, has not yet been clearly established. On the other hand, 
the phenomenon has been seen quite convincingly in lower dimensional gauge theories such 
as QED2 H, QCD3 Hi, or in QCD4 at finite temperature 0. The expected softening 
(though not flattening) of the potential due to sea-quarks has also been seen in recent work 
on three flavor QCD g. 

The reasons suggested for the failure to observe a clear signal of string breaking in zero 
temperature 4-dimensional QCD range from the inability to decouple the lowest energy 
string states at the still rather small Euclidean time extents of the measured Wilson loops 
^ to the existence of a completely new phase of the string in which breaking is completely 
invisible Q]. In the former case, it has been suggested that use of an improved string operator 
which suppresses appropriately the coupling of higher energy string states as the breaking 
point is approached is a prerequisite for exposing the desired flattening of the string tension 
. The results that we present in this paper conflrm that string breaking only appears in 
Wilson line correlators at sufficiently large Euclidean times, but demonstrate the breaking 
directly at the level of unsmeared Wilson lines and without explicitly mixing in two-meson 
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states [g[ in the initial string state. 

The three major differences between the simulations described here and previous studies 
of stringbreaking are 

1. We measure the static energy of Wilson lines in Coulomb gauge rather than the large 
area behavior of Wilson loops. The Wilson line operator in Coulomb gauge has a 
larger overlap with the lowest energy states of a static quark-antiquark source pair 
than the Wilson loop. (In the case of abelian gauge theory, the overlap for the Wilson 
line operator is perfect.) 

2. We work on physically large coarse lattices, which allows us to go to large Euclidean 
time (up to 1.2 fm) in order to project out the lowest state, which for large distance 
corresponds to a meson-antimeson pair. 

3. We work at relatively small quark mass (pion mass about 195 MeV). As emphasized 
previously |^ , the truncated determinant algorithm works perfectly well at arbitrarily 
small quark masses, as the convergence rate of the Lanczos algorithm used to extract 
the low eigenvalues depends only on the density of the infrared spectrum and does 
not dramatically deteriorate as we lower the quark mass. 

The dynamical fermion algorithm used in this paper has been described in considerable 
detail elsewhere so we will merely summarize the basic features. The hermitian quark 
Dirac operator has a completely gauge invariant spectrum which can therefore be gauge- 
invariantly split into a set of infrared modes (up to some momentum cutoff large enough 
to encompass the desired infrared virtual quark physics), while the ultraviolet modes are 
very accurately modelled by a local gauge-invariant pure-gauge action |^ which for large 
distance physics mainly results in a coupling rcnormalization. In the truncated determinant 
algorithm (TDA), we perform a simulation of the theory including exactly all the infrared 
eigenvalues in the quark determinant up to the chosen cutoff (for a detailed description of 
this approach see |^). The Lanczos technique used to extract the low eigenvalues does not 
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suffer from rapidly increasing convergence time even at very low quark masses (in contrast 
to HMC algorithms, where the quark inversions become prohibitively expensive in this 
limit), so we are able to work essentially at kappa critical. Of course, string breaking is 
expected to set in earlier for light dynamical quarks, so the TDA approach has a natural 
advantage over other dynamical QCD schemes for this problem. The other important feature 
of the calculations described in this paper is the use of large coarse lattices with an O(a^) 
improved gauge action to restore the rotational invariance of the measured static energies. 
Specifically, we give results for simulations performed on 6"* lattices with the improved action 
S =3.7[1.0(plaq) + 1.04(trt)], using the notation of [|lOj. Although there is considerable scale 
uncertainty on such a course lattice, we estimate the lattice spacing for this theory is about 
0.4 fm. We have performed the simulations at k=0.2050 corresponding to 111^^ ~ 195 MeV. 
Preliminary results from these simulations have been reported earlier [^ . 

In Section 2 we describe the computational and statistical issues underlying our results. 
Some further details of the computational load required by the TDA method are discussed. 
We describe the equilibration of our configurations in the TDA simulations, and the autocor- 
relation data underlying our error analysis are given. In Section 3 results for the correlators 
of Wilson lines at Euclidean times 0.4, 0.8, and 1.2 fm (1,2 and 3 lattice spacings) are re- 
ported. The flattening of the static energy associated with string breaking is finally visible 
for T~ 1 fm. In Section 4 we summarize our conclusions and indicate ongoing calculations. 
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2 TDA Simulations on large coarse lattices 

In the TDA approach to dynamical QCD, the quark determinant ^{A) = det{j5{I^{A) — 
m)) = det(ff) is gauge-invariantly spht into infrared and ultraviolet parts 

V{A)=Vir{A)Vuv{A) (1) 

by introducing a cutoff Acut on the absolute magnitude of the eigenvalues Aj of the hermi- 
tian Dirac operator H. These eigenvalues are gauge-invariant generalizations of the quark 
offshellness (i.e. for ^ — *■ 0, Aj — > ±i/p2 + m?)) and the cutoff is chosen to include as much 
as possible of the low energy structure of the unquenched theory while leaving the fluctua- 
tions of hiDiii{A) (which is included exactly in the Boltzmann measure of the simulation) 
of order unity after each sweep updating all gauge links. Fortunately, this choice is possible 
on lattices of large physical size as well as for light quark masses close to the critical value. 
In the simulations reported in this paper, the cutoff Acut is chosen at about 420 MeV. On 
the lattices generated this corresponds to including the lowest 840 eigenvalues of H. These 
eigenvalues are extracted by a Lanczos procedure for each trial gauge configuration gen- 
erated using an improved gauge action and the new configuration is then subjected to a 
Metropolis accept/reject step based on the change in Nfig^ylnV iii{A) (we have used A^flav=2 
flavors of degenerate light quarks in the simulations). Detailed balance in this procedure is 
ensured by using a random link choice procedure in the pure gauge update step. We have 
stored gauge configurations after ten combined gauge-update -|- metropolis determinant 
accept/reject steps. (The metropolis step had a typical acceptance rate of 50%.) 

The computational load in these simulations is completely dominated by the extraction 
of the infrared quark eigenmodes (typically to seven or eight place accuracy). For example, 
on a 6^ lattice, the gauge update takes a few seconds while the calculation of "Dm takes 
about 9 minutes on an Athlon 650- MHz processor. The string breaking results reported in 
Section 3 were performed on such a processor and required about 2.0 processor-months to 
accumulate 1000 configurations once equilibration is reached. 
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As mentioned previously, we work on coarse 6* lattices but with O(a^) improved gauge 
action. Following Alford et al jl^, we improve the gauge action with a single additional 
operator, with coefhcients tuned to optimize rotational invariance of the string tension 



SiU) = /3piaq^ jRcTr(l-C/piaq) 

plaq 

+ /3trt^^ReTr(l-C/trt) (2) 

trt 

where "trt" refers to a 8 link loop of generic structure (+x,+y,+x,-y,-x,+y,-x,-y) (the 



"twisted rectangle" of Ref|lO|). With the choices /3piaq=3.7, /3trt = l-04/9piaq, the quenched 
static quark potential becomes a smooth function of lattice radial separation even on 
these very coarse lattices, with lattice spacing a ~ 0.4 fm (for the unquenched theory). As 
we do not improve the quark action, the lattice spacing quoted here is determined by match- 
ing the initial linear rise of the string tension to a physical value. The results given below 
show that the restoration of rotational invariance survives reasonably well the introduction 
of the quark determinant, so that we have not found it necessary to retune the pure gauge 
action. 

To maximize our chances of seeing string breaking within the spatial limitations of the 
lattices being used we have chosen a kappa value corresponding to a rather light pion, 
namely k=0.2050, corresponding to a pion mass of about 195 MeV (or 0.39 in lattice units). 
This does not seem to result in a serious loss of acceptance at the level of individual Monte 
Carlo steps, but the equilibration process (starting from a quenched initial configuration) 
is definitely slower in comparison to TDA simulations performed with heavier quarks on 
lattices of smaller physical size Q . The sequence of infrared determinant values (specifically, 
hiDiji(A)) generated in a TDA simulation starting from a quenched configuration is shown 
in Fig.l. It appears that about 7000 Monte Carlo sweeps were needed to reach equilibrium 
(corresponding to about 1.5 650-MHz Athlon processor months). The results reported in the 
next section were based on the 1000 saved configurations between sweeps 7000 and 169900 
(see Figure ||) . 
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Figure 1: Equilibration of the quark determinant in TDA simulation of a coarse 6^ 
lattice 
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Figure 2: Typical autocorrelation curves of Wilson line correlators (R=l,3; T=3) 
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Given the rather long times needed for equilibration, the issue of autocorrelation of mea- 
sured quantities in the post-equilibrium configurations naturally assumes great importance. 
We have performed a careful study of the decorrelation rate of all the Wilson line correlators 
used to extract the static energy plots given in Section 3 below. Thus, if W{t) is a Wilson 
line correlator measured for the r'th saved configuration, an autocorrelation function can 
be computed as the ensemble average 



The autocorrelation time can then be read off from the exponential decay of C (r) . Typical 
autocorrelation curves are shown in Figure ^, which shows the decorrelation of Wilson line 
correlators spatially separated by 1 or 3 lattice units and of Euclidean time extent 3, mea- 
sured on an ensemble of 300 successive saved configurations. For i? = 1, an exponential fit 
gives an autocorrelation time of 7.9 (in units of 10 Monte Carlo steps), while the area under 
the autocorrelation curve gives 8.4 for the same quantity. For R — 3, the autocorrelation 
time is of order unity (1.5, from an exponential fit). In general, autocorrelation times range 
from about unity for the largest loops to on the order of 20 for the smallest. We have there- 
fore assumed that line correlators from successive bins of 20 (or more) saved configurations 
are effectively decorrelated for all relevant R and T. 

We have calculated the standard light meson and baryon spectrum, the eta prime, and 
the heavy-light mesons. We will report on these results in a future paper. The ground state 
heavy light meson mass is relevant to the static energy analysis. As shown in Figure |^, a 
consistent mass value was obtained for all the various fitting time windows. The mass fit 
using T = 1 - 5 was tuhl = 0.811 ± .007. 




(3) 
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Figure 3: Various fits to ground state heavy-light meson mass. The start time and 
end time are represented as the decimal start.end 



3 Static Energy Results 

The static energy of a quark-antiquark pair is calculated using the Euclidean time evolution 
of a color-singlet pair in Coulomb gauge. The relevant correlator is 



W{R,T) =< *(0,r)^'(i?,T)^'(i?,0)*(0,0) > 



coul 



(4) 



where ^' is an inhnitely massive quark held and the correlator is evaluated in Coulomb 
gauge. The static energy V{R) is then defined as the energy of the lowest state coupled by 
^(^,0)^'(0,0) to the vacuum, i.e. 



W{R,T) ^ Cexp{-TV{R)), T oo 



(5) 



or equivalently 



ViR)^ hm mW:^ 
T^oo W{R,T) 



(6) 
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On the lattice the correlator (4) is evaluated as an ensemble average of Tr(L(0, T)L\R, T)), 
where L{R,T) is a Wilson line- i.e. a product of T adjacent link matrices in the temporal 
direction at spatial location R. As mentioned previously, there is a distinct advantage to the 
use of Coulomb gauge Wilson line correlators over Wilson loop expectations in the study 
of stringbreaking effects. The Wilson loop necessarily involves additional contributions 
from intermediate states containing transverse gluons which are absent in the Coulomb 
gauge Wilson line correlators. In perturbation theory these states correspond to diagrams 
in which gluons are exchanged between the top and bottom horizontal (i.e. fixed time) 
portions of the Wilson loop. In quenched abelian theory, for example, the correlator (4) is a 
pure exponential W{R,T) cx e^^^^''"^, whereas the corresponding Wilson loop expectation 
is proportional to Q-i^WT+v(T)R) ij-npiyi^g the presence of excited states. It is therefore 
reasonable to expect that stringbreaking in unquenched QCD will emerge more rapidly (i.e. 
at smaller Euclidean times T) if the Coulomb gauge correlator (4) is used. 

Since the calculation is done on symmetrical 6^ lattices, we actually obtain four sets of 
correlators for each gauge configuration, obtained by successively gaugefixing to Coulomb 
gauge with the time direction chosen as each of the four original Euclidean spacetime di- 
rections. Exploiting this degeneracy provides a valuable increase in the statistics, as we 
effectively have an ensemble of 4000 correlators from the original set of 1000 gauge configu- 
rations. The autocorrelation times described in the preceding section were obtained however 
by preaveraging the four correlator sets coming from each configuration. 

In Figure ^ the static energy a,t T — 2 (circles) is compared with the similar results for 
two other theories. First, we show (down triangles) the T — 2 static energy for the quenched 
theory (with identical improved gauge action (2) but with the determinant contribution 
switched off), obtained from an ensemble of 8000 configurations, but otherwise analysed 
identically to the dynamical configurations (the untuned quenched theory). Both a shift in 
scale at short distances and the deviation at longer distances are clearly observed. Second, 
the static energy a.t T — 2 for the original Alford et. al. ||l^ unquenched action is plotted 
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Figure 4: Static energy from Wilson line correlator (T=2) for unquenched QCD and 
quenched (untuned and tuned) QCD. 
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Figure 5: Static energy from unquenched QCD Wilson line correlators (T=l,2,3). 
The calculated heavy- light meson pair threshhold is also indicated. 



(up triangles). This curve was used to tune our action to obtain an approximately equal 
initial slope. 

In Figure]^ we show the static energy curves for T=l,2 and 3 in (6) respectively. For 
comparison the heavy-light meson pair production threshhold is also shown. As we are 
working on a large coarse lattice (lattice spacing a ~0.4 fm) 3 temporal lattice spacings 
already represents a fairly large Euclidean time T ~ 1.2 fm and the signal to noise ratio 
at larger spatial distances has clearly deteriorated substantially at T=3. The shortest time 
evolution, T=l shown in Fig(3), on the other hand is still contaminated by higher energy 
states and the flattening of the string tension is certainly not visible. However, at T=3 the 
string tension appears to flatten out convincingly for distances R > 1 fm (=2. 5a). (Further 
simulations are in progress to substantially reduce the statistical errors in this regime.) 
Moreover, the smoothness of the potential curve within the statistical errors suggests that 
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the improved action terms are doing a good job of restoring rotational invariance on this 
very coarse lattice. 

The errors in the static potential are obtained most simply by computing a separate 
error on the line correlators W{R, T) and W{R, T — 1) in (6), using the measured standard 
deviation and an autocorrelation time extracted separately for every loop size R, T. The 
errors for the log ratio in (6) can then be obtained by combining the numerator and denom- 
inator errors in quadrature. This approach however almost surely yields an overestimate of 
the true errors, as W{R,T) and W{R,T — 1) tend to be positively correlated, reducing the 
variance in the ratio. Line correlators extracted from individual configurations on a small 
lattice tend to be extremely noisy (especially for the large loops of interest here) so to exam- 
ine this correlation we have performed the error analysis by binning the 1000 configurations 
into 40 sets of 25 consecutive configurations. The average line correlators in each bin are 
then essentially decorrelated and a straightforward jackknife analysis can be performed on 
the ratios in (6). This is the approach used to obtain the errors displayed in Figure |^. 
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4 Conclusions 



The analysis of dynamical two-flavor configurations on large coarse lattices obtained by 
the truncated determinant simulation method provides evidence for string breaking in zero 
temperature 4-dimensional QCD. One clear advantage of the method essential for its success 
in this case is the ability to generate a sufficient number of equilibrated and decorrelated 
dynamical configurations even for light quark mass on a physically large lattice, where 
standard HMC simulations would encounter computational problems. Evidently, as is the 
case for QED2 the infrared quark modes included in the determinant in the TDA contain 
all the essential physics of string breaking. 

Simulations on larger coarse lattices (8^ lattices with the same action (2)) as well as at 
different kappa values for the two light quark flavors are in progress. In particular, we hope 
to study in detail the hairpin amplitudes relevant to the eta-prime mass. Here also, the 
essential physics should be obtainable from amplitudes measured in a TDA simulation Q. 
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